*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*	This file generates VIV plots (Figure 1 and A1)
*	----------------------------------------------------------------------------
*	IN: 	vam_analysis_file
*
*	OUT:	//name of output
			local outname 		"${city}_VIV_\`y'"

*	----------------------------------------------------------------------------
	args sample sch_res ptype bw
	tokenize `sample', parse("_")
	local years "`5'`6'`7'`8'"

	* settings
	local bins     20
	local vams 	   OLSVAM_func OLSVAM_flag POLSVAM_func POLSVAM_flag

	* estimate test stats
	local estimate 0

	* plot individual figures
	local plot     1

	* combine graphs
	local combine  1

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	cap log close
	log using "${log}${city}_viv.smcl", replace

*	load

	if `estimate' {

	foreach y in m e {

		u "${builddata}${city}_analysisfile_vam`sample'`sch_res'`bw'.dta", clear

	*	OLS-VAM

		foreach model in func flag{
			gen OLSVAM_`model'    = OLSVAM_`y'_`model'
			gen b_OLSVAM_`model'  = b_OLSVAM_`y'_`model'
		}
		foreach model in func flag{
			gen POLSVAM_`model'   = POLSVAM_`y'_`ptype'_`model'
			gen b_POLSVAM_`model' = b_POLSVAM_`y'_`ptype'_`model'
		}

		foreach vam in `vams'{

			preserve

			//bin and set vars
			ivset, ptype(`ptype')
				local S=r(S)
			vambin if `S', n(`bins') ptype(`ptype') vam(b_OLSVAM_`y'_flag)
			ivset, y(`y') ptype(`ptype') bin
				local Z=r(Z)
				local p=r(P)
				local S=r(S)
				local Y=r(Y)

		*	forecast

			eststo `vam': ivreg2 `Y' (`vam' = `Z') `p' $iv_controls if `S', partial(`p' $iv_controls) small
			*scalar forecast_coef = _b[`vam']
			local forecast_coef = _b[`vam']

			test `vam' == 1
			local forecast = r(F)

			local omnibus = (e(N)-e(df_m)-1)/e(N) * e(j) + `forecast'
			local omnibus_df = e(jdf)+1
			local omnibus_p = chi2tail(`omnibus_df',`omnibus')


		*	omnibus over-id

			scalar overid = (e(N)-e(df_m)-1)/e(N) * e(j)
			scalar omnibus = overid + `forecast'
			scalar omnibus_p = chi2tail(e(jdf)+1,omnibus)

		***	first-step estimates

			mat toplot = J(`=abs(`bins')',8,.)

		*	reduced form
			keep if `S'

			//estimate (cluster in suest later)
			qui reg `Y' `Z' `p' $iv_controls if `S'
			mat b = e(b)
			mat b = b[1,1..`bins']
			mat T = I(`bins') - J(`bins',`bins',1/`bins')
			mat b_raw = b
			mat b = b*T

			//store
			mat toplot[1,1] = b'
			mat toplot[1,6] = b_raw'

		*	first stage

			//estimate (cluster in suest later)
			qui reg `vam' `Z' `p' $iv_controls if `S'
			mat b = e(b)
			mat b = b[1,1..`bins']
			mat b_raw = b
			mat b = b*T

			//store
			mat toplot[1,2] = b'
			mat toplot[1,7] = b_raw'

		*	testing for equality (rf = fs)

			qui sureg (`Y' `Z' `p' $iv_controls) (`vam' `Z' `p' $iv_controls)
			local i 0
			forval bin=1/20 {
				qui test [`Y']_b[binoffer_`bin'] = [`vam']_b[binoffer_`bin']
				mat toplot[`++i',8] = r(p)
			}

			//add school code to matrix (for merge)
			if `bins'==0{
				local i 0
				foreach offer of varlist offer_*{
					local sch = substr("`offer'",7,.)
					mat toplot[`++i',4] = `sch'
				}
			}

			//bin number
			if `bins'!=0{
				local i 0
				forval j =1/`bins'{
					mat toplot[`++i',5] = abs(`bins')-`j'+1
				}
			}

			//load in matrix of first-steps
			clear
			svmat toplot
			rename toplot1 rf
			rename toplot2 fs
			rename toplot3 t
			rename toplot4 sch
			rename toplot5 n
			rename toplot6 rf_raw
			rename toplot7 fs_raw
			rename toplot8 p
			*tostring n, replace

			//store forecast and overid
			g forecast_coef = `forecast_coef'
			g omnibus_p		= `omnibus_p'

			* save intermediate file
			sa "${builddata}${city}_`vam'_`y'_`ptype'_viv.dta", replace

			restore

			}
		}
	}


	if `plot' {

		local vams POLSVAM_flag_m OLSVAM_func_m OLSVAM_flag_m POLSVAM_func_m POLSVAM_flag_m ///
				   OLSVAM_func_e OLSVAM_flag_e POLSVAM_func_e POLSVAM_flag_e

		foreach vam in `vams' {

		* load
		u "${builddata}${city}_`vam'_`ptype'_viv.dta", clear

		* demean fs and scale rf by same factor
		qui su fs_raw
		replace fs = fs_raw - r(mean)
		replace rf = rf_raw - r(mean)

		***	plot

			//parameters
			if "${city}" == "DEN" {
				local texty_fc = -0.34
				local texty_om  = -0.39
				local textx = 0.03
				local min = -0.4
				local max =  0.4
				local gap = 0.2
			}
			if "${city}" == "NYC" & !inlist("`vam'","OLSVAM_func_m","OLSVAM_func_e") {
				local texty_fc = -0.34
				local texty_om  = -0.39
				local textx = 0.03
				local min = -0.4
				local max =  0.4
				local gap = 0.2
			}
			if "${city}" == "NYC" & inlist("`vam'","OLSVAM_func_m","OLSVAM_func_e") {
				local texty_fc = -0.31
				local texty_om  = -0.39
				local textx = 0.235
				local min = -0.4
				local max =  0.8
				local gap = 0.2
			}
			if "${city}" == "NYCms" & !inlist("`vam'","OLSVAM_func_m","OLSVAM_func_e") {
				local texty_fc = -0.34
				local texty_om  = -0.39
				local textx = 0.03
				local min = -0.4
				local max =  0.4
				local gap = 0.2
			}
			if "${city}" == "NYCms" & inlist("`vam'","OLSVAM_func_m","OLSVAM_func_e") {
				local texty_fc = -0.31
				local texty_om  = -0.39
				local textx = 0.13
				local min = -0.4
				local max =  0.6
				local gap = 0.2
			}
			local window rf>`min' & rf<`max' & fs>`min' & fs<`max'
			local xtitle "First stage effect on OLS VAM"
			local ytitle "Reduced form effect on test scores"
			if strpos("`vam'", "OLSVAM_func") >= 1 local title "Uncontrolled"
			if strpos("`vam'", "OLSVAM_flag") >= 1 local title "Conventional"
			if strpos("`vam'", "POLSVAM_func") >= 1 local title "Risk only"
			if strpos("`vam'", "POLSVAM_flag") >= 1 local title "RC VAM"

			//added text
			qui su forecast_coef
			scalar forecast_coef = r(mean)
			local forecast_coef = r(mean)
			qui su omnibus_p
			scalar omnibus_p = r(mean)

			local forecast_text "Forecast coef.: `:di %4.3f forecast_coef'"
			if `forecast_coef' > 1 local forecast_text "Forecast coef.: `:di %4.2f forecast_coef'"
			local omnibus_text  "Omnibus p-val.: `:di %4.3f omnibus_p'"

			//coordinates for plotted lines
			gen eqline = `min' in 1
			replace eqline = `max' in 2
			if forecast_coef < 1{
				gen slopex = eqline
				qui su rf
				gen slopey = slopex*forecast_coef + r(mean)
			}
			else{
				gen slopex = `min'/forecast_coef in 1
				replace slopex = `max'/forecast_coef in 2
				qui su rf
				gen slopey  = slopex*forecast_coef + r(mean)
			}
			local eqline line eqline eqline, lcolor(gs9) lwidth(thin) lpattern(dash)
			local fcline line slopey slopex, lcolor(gs4) lwidth(medthin) lpattern(solid)

			* bound fit line so plot is square
			qui reg rf fs
			local intercept = _b[_cons]
			local minrange = (`min' - `intercept') / _b[fs]
			local maxrange = (`max' - `intercept') / _b[fs]
			if `minrange' < `min' local minrange = `min'
			if `maxrange' > `max' local maxrange = `max'

			local legendlabs label(3 "45 degree line") label(4 "Forecast coef. line")
			local mlabopts mlabcolor(gs4) mlabp(4) mlabsize(vsmall) mlabgap(.4)

			*	draw
			twoway 	scatter rf fs if p<=.10 & `window', msymbol(d) color(gs4) `mlabopts' || ///
					scatter rf fs if p>.10 & `window', msymbol(dh) color(gs4) `mlabopts' || ///
					`eqline' || `fcline' ///
					subtitle("`title'", size(large)) xtitle("`xtitle'", color(gs3) size(large)) ytitle("`ytitle'", color(gs3) height(10) size(large)) ///
					xlab(`min'(`gap')`max', nogrid) ylab(`min'(`gap')`max', nogrid) ///
					text("`texty_fc'" "`textx'" "`forecast_text'", placement(right) size(*1.1) color(gs6)) text("`texty_om'" "`textx'" "`omnibus_text'", placement(right) size(*1.1) color(gs6))	///
					legend(order(4 3) `legendlabs' rows(1) position(6) span region(lc(white)) size(small) bmargin(small)) ///
					yline(0, lcolor(gs15)) xline(0, lcolor(gs15)) ///
					aspectratio(1) scale(.95) graphregion(color(white)) scheme(s2mono) ///
					name(${city}_`vam'_`ptype', replace) nodraw

			graph save "${figures}${city}_`vam'_`ptype'.gph", replace

		}
	}


	if `combine' {

		foreach y in m e {

		* combine
		grc1leg ${city}_OLSVAM_func_`y'_`ptype' ${city}_OLSVAM_flag_`y'_`ptype' ${city}_POLSVAM_func_`y'_`ptype' ${city}_POLSVAM_flag_`y'_`ptype', ///
			altshrink rows(1) colfirst graphregion(color(white)) name(combined, replace) legendfrom(${city}_OLSVAM_func_`y'_`ptype') plotr(m(small))
		graph display combined, ysize(2) xsize(4)

		graph export "${figures}`outname'.pdf", replace

		}
	}

	log close
